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Abstract. After introducing the fundamental properties of self-gravitating systems, we 
present an application of Tsallis' generalized entropy to the analysis of their thermody- 
namic nature. By extremizing the Tsallis entropy we obtain an equation of state known 
as the stellar polytrope. For a self-gravitating stellar system confined within a perfectly 
reflecting wall, we discuss the thermodynamic instability caused by its negative specific 
heat. The role of the extremum as a quasi-equilibrium is also demonstrated from the 
results of iV-body simulations. 
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1 Introduction 

In any subject of astrophysics and cosmology, many-body gravitating systems play an essential role. 
Globular clusters and elliptical galaxies, which are recognized as self-gravitating stellar systems, are 
typical examples | ll'2mp4 |. Several beautiful works on the thermodynamics of self-gravitating systems 
5 , 6 have shown peculiar features such as a negative specific heat and the absence of global entropy 
maxima, which is referred to as the gravothermal catastrophe. (For a general introduction to this 
subject, see 0.) Furthermore, the long-range nature of the gravitational interaction makes a discussion 
of the relationship between non-extensive thermostatistics the self-gravitating systems tempting. 

Recently, a new framework for thermodynamics based on Tsallis' non-extensive entropy was pro- 
posed 8 . It has been applied extensively to deal with a variety of interesting problems to which stan- 
dard Boltzmann-Gibbs statistical mechanics can not be applied |§lll(J| . The study of self-gravitating 
stellar systems has been one of the most interesting applications of Tsallis' framework of thermostatis- 
tics (see, e.g., |11II12II13| 1. Here, some of the progress in its application to stellar systems )14II15II16II17| 
will be reviewed. Although its dynamics is complicated in general, if we impose spherical symmetry 
on the system, its treatment is considerably simplified due to the well-known 1/r 2 behavior of the 
gravitational force. Furthermore the spherical system still keeps the remarkable nature of a negative 
specific heat Thus, self-gravitating stellar systems seem to be a desirable testing ground for 

Tsallis' non-extensive thermostatistics. 

This paper is organized as follows. In section 2, we briefly review the principal characteristics of 
self-gravitating systems. In particular, the standard treatment for the gravothermal catastrophe based 
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on the Boltzmann-Gibbs entropy is explained. Then, its extension to the Tsallis entropy is discussed in 
section 3. The properties of states of the stellar system that extremize the Tsallis entropy are clarified. 
In section 4, the interesting role of the above states as quasi-equilibria is discussed from the results of 
numerical simulations. Finally, section 5 is devoted to discussion and conclusions. 



2 Some basic properties of the self-gravitating systems 

In order to get an intuitive understanding of the evolution of many-body self-gravitating systems, let 
us consider a very simple situation: the circular motion of a particle of mass m and velocity v at a 
radius r in a spherical mass distribution with a constant density po: 

v 2 GmM(r) 4?r „ , , 

m— = j±±, Mr=— r 3 p , (2.1 

r r z 6 

where G is the gravitational constant and M(r) is the mass contained within the sphere of radius r. 
By means of the orbital period T, we define the dynamical time of the system (see, e.g., p. 57 of pQ) 
as 

tdya = i = V^fe' (2 ' 2) 

which is recognized as the characteristic time scale of more general self-gravitating systems with mean 
density pq. A distribution of the self-gravitating system evolves with the time scale H2.2L so this 
quantity is called the dynamical time. As explained in the Appendix, the many-body self-gravitating 
system has another timescale, i.e. the relaxation time ^S] (see also chapter 4 of pQ), 

0.1N 

~ TniV <dyn 

where N is the number of particles in the system. This time scale represents the relaxation processes 
due to scattering by the gravitational interaction between each pair of particles. 

In a usual system, we believe that a distribution of particles approaches the so-called thermal 
equilibrium state within the relaxation time l|2.3[l . However, this is not necessarily the case for self- 
gravitating systems, since they have the peculiar property of a negative specific heat, which is explained 
through the following discussion. For the circular orbit obeying (|2.1|) . the virial theorem is easily derived 
as 

2K + U = 0, (2.4) 

where K = mv 2 /2 is the kinetic and U = —GmM(r)/r the potential energy. In general, for self- 
gravitating many-body systems consisting of N particles, a virial relation identical to (|2.4|) holds 
between long-time-averaged values of the total kinetic and the total potential energy (see, e.g., chapter 
8.1 of 1 ). Thus, the total energy E can be expressed as 

E = K + U = -K. (2.5) 

If we introduce a temperature to characterize the total kinetic energy as K = 3NkT/2, where k is the 
Boltzmann constant, then l|2.5[) tells us that the specific heat of a self-gravitating system is negative, 
since 

„ dE 3Nk n 

C= — = < 0. (2.6 

dT 2 K ' 

which suggests the existence of a thermodynamic instability proceeding with the relaxation timescale 
(l2~31) . 

For a more precise discussion of the thermodynamic properties of self-gravitating systems, Antonov 
[5] and Lynden-Bell and Wood [Sj considered a system confined within a spherical adiabatic wall of 
radius r e . This means that the N particles in this system interact via Newtonian gravity and bounce 
elastically from the wall. In order to investigate the equilibrium states of the system, the Boltzmann- 
Gibbs entropy of the system, 

3bg = - ff In / d 6 r, (2.7) 
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Fig. 2.1. Energy-radius-density contrast (left) and radius-temperature-density contrast (right) relationships for the 
isothermal stellar system 



is introduced, where f(x, v) is the distribution function in phase-space. Here, the phase space measure 
is defined as 

K d 3 xd 3 v , , , 
d 6 r = — , h 3 = (l v ) 3 (2.8) 

where h 3 denotes the phase-space element with unit length Iq and unit velocity vq. Hereafter, we set 
the Boltzmann constant to be unity, i.e. k = 1. Under the condition that the total mass M and the 
energy E of the system are kept constant, we maximize the entropy Sbg : 



SS BG ~ aSM - (38 E = 0, 



(2.9) 



where a and fj are Lagrange multipliers. Through this procedure, the equation of state of the system 
turns to be isothermal [S]: 



, 1 9 d 3 v 



p(x) 



f 



d 3 v 



P( x ) 



(2.10) 



where p and p are the pressure and density, respectively The thermal equilibrium states of ordinary sys- 
tems are homogeneous; i.e. they have uniform density. However, this is not the case for self-gravitating 
systems. For the spherical systems we are discussing here, the magnitude of the density and pressure 
should increase as the radius r decreases in order to support the system against its self-gravity It is 
shown that a series of equilibrium states is well parameterized by the value of the density contrast 
D = Pc/pe, where p c is the central density and p e = p{r e ) is the density at the wall. In the next 
section, we will present the procedures for calculating the equilibrium state in more general situations. 
In Fig. 12.11 (left), we plot the dimensionless quantity A = —r e E/GM 2 as a function of the density 
contrast D. Since the solid curve in this figure represents states which extremize the entropy, we note 
that no equilibrium state is attained when the radius of the wall r e is larger than the critical radius 

Gil/ 2 

r c = 0.335--^y, (2.11) 

which corresponds to the state with the critical density contrast D c = 709. Furthermore, it can be 
shown, from a turning-point analysis for a linear series of equilibria (6lll9ll20ll21| . that, along the curve 
\(D) derived from the condition SSbg = 0, all states with D > D c are unstable. Also, the explicit 
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evaluation of the eigenmodes for the second variation of the entropy S 2 Sbg leads to the same results 
[7lll9| . We call this phenomenon, i.e. the absence of thermal equilibrium in self-gravitating systems 
when r > r c and/or D > D c , the gravothermal catastrophe. 

Heuristically, this instability is explained by the presence of the negative specific heat as follows. 
In a fully relaxed gravitating system with a sufficiently larger radius, the negative specific heat arises 
at the inner part of the system, where we have Cv,inncr < 0, while the specific heat at the outer part 
remains positive (Cy, outer > 0), since one can safely neglect the effect of self-gravity. In this situation, 
if a tiny heat flow is momentarily supplied from the inner to the outer part, both the inner and outer 
parts get hotter after the readjustment of the system. Now imagine the case Cy, outer > I Cy, inner I- The 
outer part has so much thermal inertia that it cannot heat up as fast as the inner part, and thereby 
the temperature difference between inner and outer parts increases. As a consequence, the heat flow 
never stops, leading to a catastrophic temperature growth. 

Fig. l2.ir right'l plots the dimensionless inverse temperature r\ = GMf3/r e with respect to the density 
contrast. We can evaluate the specific heat at constant volume, 




(2-12) 



as shown in Fig l2.2l 

In contrast to the naive discussion for i|2.6[l based on the virial theorem, the actual specific heat 
(I2.12|) changes sign many times. We observe that a self-gravitating system confined by a wall of smaller 
radius has a positive specific heat. The reason is as follows. It can be shown that the virial theorem 
for a system surrounded by a boundary leads not to (|2.4|) but to 

IK + U = 3P e V, (2.13) 

where P e is the pressure at the wall and V the volume of the system (see, e.g., eq.(29) of QI]). For 
the case of a smaller radius, the surface energy dominates the gravitational energy of the bulk so that 
the system behaves as a normal one; i.e. Cy > 0. As the radius r e increases, the gravitational energy 
becomes significant. Eventually, the specific heat changes its sign, Cy < , when the value of the 
density contrast becomes D = 32.2, which corresponds to a peak of the rj(D) curve (the left side of 

FigEH)- 

At the onset of the gravothermal catastrophe, D c = 709, the specific heat vanishes, and Cy > in 
the unstable segment of the X(D) curve, D > 709. This fact implicitly indicates a role of the negative 
specific heat in existence of the thermodynamic instability in self-gravitating systems. Following the 
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previous heuristic argument, above the critical point (D > D c = 709), the specific heat of the total 
system becomes positive (Cy = Cy, inner + Cy, outer > 0) equivalently, 

Cy.outer > — Cy.inner > 0, (2-14) 

which means the outer normal part of the system has a larger heat capacity, so that it cannot catch 
up with the temperature increase of the inner part. 

So far, we have discussed some interesting features of self-gravitating systems. Among these prop- 
erties, the peculiarity of the system arising from the negative specific heat has been studied in terms 
of the Boltzmann-Gibbs entropy (|2.7J) . In the next section, this sort of the analysis will be extended to 
that based on the non-additive entropy, i.e. the Tsallis entropy. However, before closing this section, 
we briefly summarize the subsequent investigation of the gravothermal catastrophe. 

In this paper we adopt the thermodynamic approach in terms of the entropy and/or the free energy 
in order to investigate instabilities of self-gravitating systems. Although this method has the advantage 
of simplicity, we cannot deal with dynamical aspects of the instability appearing within the relaxation 
timescale l|2.3l) . Typically, each particle (star) in a self-gravitating system has a long mean free path 
compared with its system size. In other words, during its travel across the system, the the motion of a 
particle is driven by the gravitational potential of the whole system and the particle experiences quite 
a few two-body encounters, so that the relaxation time t Te \ Ij2.3|l is much longer than the dynamical one 
t&y a (E3- This means that we cannot rely on the usual notion of local equilibrium, which is realized 
in cases with a short mean free path. Thus, the Fokker- Planck equation (see, e.g., chapter 8.3 of £Q), 
which can be derived from the expansion of a collision integral in the Boltzmann equation with respect 
to the momentum transfer due to a two-body collision, is a useful approach for studying the long-term 
evolution of a self-gravitating system. By means of the linearized Fokker-Planck equation, eigenvalues 
and eigenfunctions which describe the time evolution of perturbations around the thermal equilibrium 
have been obtained |22| . 

The gas model in which a star (particle) in the self-gravitating many-body system is replaced by 
a gas particle, is an alternative approach. Although the gas model loses some of the basic properties 
of a the self-gravitating many-body system, we note that it still shows the interesting aspects of 
the gravothermal catastrophe. In this model, the local relaxation time is shorter than the dynamical 
timescale. Thus, we can introduce the usual notions of local equilibrium and the local temperature, 
which make our analysis much easier than in the stellar system. The local specific heat is evaluated to 
show that the condition (|2.14|) characterizes the onset of the gravothermal catastrophe [23] . The time 
evolution of a linear perturbation around the equilibrium gas with heat transfer has been investigated 
to show that the proposed mechanism for the instability based on the negative specific heat actually 
works in self-gravitating gaseous systems |24| . 



3 Analysis of gravothermal catastrophe based on Tsallis' generalized entropy 

3. 1 Tsallis entropy and stellar polytrope 

Here we shall discuss a generalization of the analysis in previous section based on the Boltzmann-Gibbs 
entropy l|2.7|l to the case of the Tsallis entropy. The details are given in the literature |14U15lll6j . The 
the energy and mass of a one- particle distribution f(x,v) are respectively expressed as 

E = K + U = y ji« 2 + i <£(*)} f(x,v)d 6 T, (3.1) 

M = j f(x,v) d 6 r, (3.2) 
where the quantity <P(x) is the gravitational potential given by 

^' = - G /|^^' < 3 - 3 > 

As for generalization, the most crucial problem is the choice of the statistical average in non- 
extensive thermostatistics. Following a seminal paper [llj . the analyses in a couple of our papers |14l 
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115) have been done by utilizing the old Tsallis formalism with the standard linear mean values (see also 
25,26 for comparative works). On the other hand, a more sophisticated framework by means of the 
normalized g-expectation values has recently been presented (27028| . As several authors have advocated, 
the analysis using normalized g-expectation values is thought to be essential, since the undesirable 
divergences in some physical systems can be eliminated safely when introducing the normalized q- 
expectation values. Furthermore, non-uniqueness of the Boltzmann-Gibbs theory has been shown by 
using normalized g-expectation values 1291 . In this paper, we mainly report our investigation based 
on normalized g-expectation values. However, this does not imply that all the analyses with standard 
linear means or un-normalized ^-expectation values lose the physical significance. Actually, we observed 
some differences between results based on two frameworks, which will be mentioned at the end of this 
section. 

In the new framework of Tsallis' non-extensive thermostatistics, all the macroscopic observables of 
the quasi-equilibrium system can be characterized by the escort distribution, but the escort distribution 
itself is not thought to be fundamental. Rather, there exists a more fundamental probability function 
p(x, v) that quantifies the phase-space structure. With the help of this function, the escort distribution 
is defined and the macroscopic observables are expressed as normalized ^-expectation value as follows 
(see, e.g., E3EBJ): 

escort distribution : P q (x, v) = — - , (3.4) 

J d 6 T{ P (x,v)} q 

normalized q- value : (Oi) q — J d 6 r Oi P q (x,v) . (3-5) 

Based on the fundamental probability p(x,v) 7 the Tsallis entropy is given by 

S q = -^-j J d 6 r {{p(x, v)} q - p(x, v)} . (3.6) 

Note that the probability p(x, v) satisfies the normalization condition: 

J d 6 r p{x,v) = 1. (3.7) 

To apply the above Tsallis formalism to the present problem without changing the definition of 
energy and mass (|3.1(l and l|3.2|l . we identify the one-particle distribution with the escort distribution 
P q , not the probability function p(x,v): 

f{x, v) = M M^H! ; Nq = J <? T {p{Xi v) y (3.8) 

so as to satisfy mass conservation (|3.2|l . The definitions of the Tsallis entropy (|3.6|l and the escort 
distribution 13.8fl clearly show that the Boltzmann-Gibbs entropy l|2.7|l is recovered in the limit q — * 1. 
Now, adopting the relation l|3.8l) . let us seek the extremum-entropy state under the constraints l|3.1() 
and 13.71) . The variational problem is 



■V " \ I <I"ti>- M / <P-r[\v 2 + \$)f -E 



0, (3.9) 



where the variables a and (5 denote the Lagrange multipliers. The variation with respect to the prob- 
ability p(x,v) gives the so-called stellar polytrope as the extremum-entropy state ( see a l so P-223 
offfl): 



where we define the constants A and <?o as 



0o - \v 2 - $(x) 



9/(1-9) 

(3.10) 



M f q(l-q) PMy'^ N q 
A ~W q \a(l- 9 ) + l AT/ ' <P °-/3M ( l-, ) +<£> ' (3 - H) 
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with the quantity < s > being 



J_ fd e r(^v 2 • </>; C-!.12i 



By means of the distribution function (|3.10() , we evaluate the density p{r) and the isotropic pressure 
P(r) at the radius r = \x\ as 

d 3 v 



P(r) = J -j^r f( x > v )> 



AV2nB (1^) ^ [<Po-Hr)] 1/{1 - q)+1/2 (3.13) 



and 



P(r) = J ^ j^ 2 

= ^ B (I' 1^) P [ ^ - %)] 1/(1 - 9)+3/2 - (3-14) 
Here, the function B(a, b) denotes the beta function. These two equations lead to the polytropic relation 

P(r) = K nP 1+1 ' n {r), (3.15) 



with the polytrope index n related to the parameter q by 

1 1 



1 -q 2 



(3.16) 



The explicit form of the dimensional constant K n is given in |16| . We note that the isothermal equation 
of state ()2.1Q(I discussed in the previous section arises in the limit n — > oo (q — > 1) as expected. 
In terms of K n and p, the one-particle distribution can be rewritten as 

f(x,v) = 



4V2n 5(3/2, n- 1/2) {(n + 1) K n pV"}W 

1 (3.17) 



(n + 1) if n p 1 /^ 

which agrees with the result based on the old Tsallis formalism. That is to say, the equilibrium distri- 
bution H3.17JI turns out to be invariant irrespective of the choice of the statistical averages. 

We note that the equilibrium distribution l|3.10|l with l|3.11f) contains the new quantities N q and 
< e >, which implicitly depend on the distribution function itself. In marked contrast to the result 
in the old Tsallis formalism, this fact gives rise to non-trivial thermodynamic relations as follows. By 
using the definitions of density and pressure l|3.13() and (|3.14|) . the quantity < e > becomes 

< £ > = 77 It: f d 3 x P{x) + [ d 3 x p(x) 



d 3 xP(x) - / d 3 x p(x) [<P Q - <P(x)} \ + <P 



M {2 

= J_(i 

M \ 2 

Furthermore, by using the relation <Pq — <P(x) — (n + l)(P/p) from H3.13(l and H3.14|l and substituting 
the equation . 1 1|) into the above expression, the quantity < e > is canceled and the equation reduces 
to 



d 3 xP(x). (3.18) 



As for N q , the normalization condition l|3.7(l gives us its actual expression |lfi) . 
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3.2 Gravothermal catastrophe in the stellar polytropic system 

We will specifically focus on the spherically symmetric case with the polytrope index n > 3/2 (q > 0), 
in which the equilibrium distribution is at least dynamically stable (see chapter 5 of In this case, 
the stellar equilibrium distribution can be characterized by the so-called Emden solutions (see, e.g., 
30 3^) and all the physical quantities are expressed in terms of the homology invariant variables 
(u, v), which are subsequently used in later analysis. 

We note that the one-particle distribution function l|3.10|) does not yet completely specify the 
equilibrium configuration, due to the presence of the gravitational potential, which implicitly depends 
on the distribution function itself. Hence, we need to specify the gravitational potential or density 
profile. From the gravitational potential (|3.3() . we obtain the Poisson equation 

1 d f 2 d<P(r) 



r i dry ^) = ^G P (r). (3.19) 

By combining the above equation with H3.13f) . we obtain the ordinary differential equation for <P. 
Alternatively, a set of equations which represent the hydrostatic equilibrium are derived by using 
(l3~T51) . and jPHh 

dP(r) Gm(r) 



dr 
dm(r) 
dr 



p(r), (3.20) 
i— 

4np(r)r 2 . (3.21) 



The quantity m(r) denotes the mass evaluated at the radius r inside the wall. Denoting the central 
density and pressure by p c and P c , we then introduce the dimcnsionless quantities: 

P = Pc[0(O}\ r = l ^±^ V 2 e, (3-22) 



AttG P 2 c 

which yields the ordinary differential equation 

0" + j6' + 0" = 0, (3.23) 

where the prime denotes the derivative with respect to £. The quantities p c and P c in l|3.22|l are the 
density and the pressure at r = 0, respectively. To obtain the physically relevant solution of 13.23fl . we 
use the boundary condition 

6»(0) = 1, 61' (0) = 0. (3.24) 

A family of solutions satisfying l|3.24[) is referred to as the Emden solution, which is well-known in the 
subject of stellar structure (see, e.g., chapter IV of ref.|30|\ To characterize the equilibrium properties 
of Emden solutions, it is convenient to introduce the following set of variables, referred to as homology 
invariants |3771l3"T] : 

d\nm(r) 47rr 3 p(r) £0™ . 

u =— n = n — = — aT> ( 3 - 25 ) 

ctlnr m[r) 

amr "(r) r 

which reduce the degree of l|3.23|l from two to one. Recall that we are investigating the thermodynamic 
properties of a self-gravitating system within a wall of radius r e . We can evaluate the total energy of 
the confined stellar system in terms of the pressure P e , the density p e at the boundary r e , and the 
total mass M: 

3 f Tc 2 f r " Gm(r) dm 



3 f r " r° 
E = K + U = - dr 47rr 2 P(r) - dt 
2 Jo Jq 



1 



n — 5 



r dr 
3 f GM 2 MP e 



5V 



-(n + 1)^^ } + {n-2)4irr 3 e P ( 
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Fig. 3.1. Energy-radius- 
density contrast relationship 
for the stellar polytropes with 
various indices. 



by which the dimcnsionlcss quantity A can be expressed as a function of the homology invariants at 
the wall IT1I15II16I : 



A 



Er e 
GAP 



1 



i 



+ (n- 



2) — 



(3.27) 



As was already shown in section [5] for the isothermal case, the above dimensionless quantity A has 
an important role for the analysis of the gravothermal catastrophe. Figure l3~Tl shows A as a function 
of the ratio of the central density to that at the boundary, p c /p e - We notice that the A-curves are 
bounded from above and have peaks in the case of n > 5 (right panel). On the other hand, the curves 
for n < 5 monotonically increase (left panel). It follows that the stellar polytrope within an adiabatic 
wall exhibits the gravothermal instability in the case of a polytropic index n > 5. Similarly to the 
isothermal case, the evaluation of eigenvalues for the second variation of the Tsallis entropy S 2 S q gives 
the same results as the above turning point analysis in terms of the A-curve |14U15lllf)| . 



3.3 Thermodynamic instability arising from the negative specific heat 

As discussed in the isothermal case in the previous section, thermodynamic instability in stellar systems 
is intimately related to the presence of a negative specific heat |fill32| . The evaluation of the specific 
heat is thus necessary for clarifying the thermodynamic properties. In this regard, the identification 
of the temperature in the stellar system in the polytropic case is the most essential task. 

In the new framework of Tsallis' non-extensive thermostatistics, the physically plausible thermo- 
dynamic temperature Xphys can be defined from the zero-th law of thermodynamics 33,34, 333 ■ The 
thermodynamic temperature in a non-extensive system differs from the usual one, i.e., the inverse of 
the Lagrange multiplier (3. The pseudo-additivity of the Tsallis entropy and the transitivity of thermal 
equilibria suggests 

T phys = [l + (l-q)S q ] /T 1 . (3.28) 

Here, the relation H3.28J) should be carefully applied to the present case, since the verification of zero- 
th thermodynamic law is very difficult in a stellar equilibrium system with long-range interactions. 
Furthermore, even using the new formalism, the energy E remains non-extensive due to the self- 
referential form of the potential energy (see H3.1(l and (|3.3(l l. However, the consistency of the physical 
temperature given by (|3.28(1 can be shown by an alternative argument |15Ulfi| for the self-gravitating 
stellar system, where the modified Clausius relation |33U34U35| 

dS q = -l—{l + (l-q)S q }d'Q. 
Jphys 

plays a significant role. 
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From the expression for the normalization factor N q of the escort distribution i|3.8[) . the Tsallis 
entropy l|3.t)|) of the extremum state is given by 

S q = j^- q {N q -l). 

Let us substitute this expression into (|3.28|1 . Then, the relation (|3.18|1 for N q /(3 gives us the expres- 
sion of the thermodynamic temperature T p h ys as an integral of a physical quantity, i.e. the pressure 
distribution: 

Tphys = y = J d 3 xP(x). (3.29) 

We can evaluate the above integral to represent the thermodynamic temperature in terms of the 
pressure P e , the density p e at the boundary and the total mass M as 

T phys = X — IteriPe -(n + l)^ + —\ (3.30) 

n - 5 [ p e r e J 

from which the dimensionless inverse temperature is denoted by 

GM 2 (n-5)v e , 
V = = ,\ ' 3.31 

as a function of the homology invariant at the boundary ^^j. From (|3.27|) and l|3.31|l . the specific heat 
at constant volume C v is given by 

/d\ 

dE \ 2 \d{ 



Cv= -pf^- =V 2 y^- (3.32) 



In order to compare with the behavior of the specific heat in the isothermal case of Fig. 12.21 wc 
plot C v with respect to the density contrast for several values of the polytropic index in Fig. 3.2. 

As already denoted at the end of subsection l3.2l the stellar polytropic system exhibits the gravother- 
mal catastrophe when the polytropic index n is larger than the critical value, i.e. n > 5. In the isother- 
mal case, the onset of the thermodynamic instability is characterized by the condition l|2.14|l . Figure 
13.31 shows that the thermodynamic argument based on the negative specific heat for the gravother- 
mal catastrophe can be extended to the case of the stellar polytrope, i.e. the equilibrium of the 
self-gravitating system described by the extremum of the Tsallis entropy. 

In Fig. 13.31 we depict the density profiles with respect to the dimensionless radius £ (|3.22|) . Clearly, 
profiles with index n < 5 rapidly fall off and they abruptly terminate at finite radius (left panel), while 
the n > 5 cases infinitely continue to extend over the outer radius (right panel). It follows that stellar 
polytropic systems with index n > 5 are able to possess a sufficient amount of the outer normal part 
to realize the gravothermal catastrophe. 
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Before closing this section, we mention differences between the old and new Tsallis formalisms. 
First, the relationship between the polytrope index n and theTsallis parameter q in 1(3.16(1 differs from 
the one obtained previously, but is related to it through the duality transformation, q <-> 1/q (see (14) 
in Jl] or (12) in ^H])- This property was first addressed in (^[ in a more general context, together with 
changes to the Lagrangian multiplier (3. The duality relation implies that all of the thermodynamic 
properties in the new formalism can also be translated into those obtained in the old formalism. 
Secondly , in previous studies |15| based on the standard linear means, the radius-mass-temperature 
relation and the specific heat seriously depend on the dimensional parameter for the phase element h 
1)2.8(1 . By contrast, in the present case |16| . the radius-mass-temperature relation was derived from the 
non-trivial relation 1(3.18(1 . in which no such /i-dependence appears. The resultant specific heat is also 
free from this dependence, which is a natural outcome of the new framework using the normalized q- 
expectation values. Therefore, it seems likely that the new formalism provides a better characterization 
for non-extensive quasi-equilibrium systems. 



4 Reality of the stellar polytrope as a quasi-equilibrium state 

So far, we have discussed the extension of the thermodynamic analysis to the stellar polytrope which 
is obtained by applying the variational procedure to the Tsallis entropy for the stellar self-gravitating 
system. It has been shown that thermodynamic quantities such as the specific heat are useful for study- 
ing the emergence of the gravothermal catastrophe in the stellar polytrope. In Fig. 14. l| we summarize 
our results. The stellar polytrope confined by an adiabatic wall is shown to be thcrmodynamically 
stable when the polytrope index n < 5. In other words, if n > 5, a stable equilibrium state ceases 
to exist for a sufficiently large density contrast D > -D C rit, where the critical value -D C rit given by a 
function of n is determined from the second variation of the entropy around the extremum state of 
the Tsallis entropy, 6 2 S q = [T^lfT?>j . The dotted line in Fig l4.ll represents the critical value D c rit for 
each polytrope index, which indicates that the stellar polytrope at low density contrast D < D cr - lt is 
expected to remain stable. 

The above arguments indicate that, similar to the isothermal state, the stellar poly tropic distribu- 
tion can also be regarded as an equilibrium state, since it is described by the extremal state of the 
Tsallis entropy. However, the one-particle distribution function of the stellar polytrope ((3.17(1 clearly 
shows that the velocity dispersion, 

a[r)(X W) J if y2f ^ v ^ ( 4J ) 

depends on the radius r. Only in the isothermal case, n — > oo, is a spatially constant. Thus, it 
is expected that a gradient of the velocity dispersion is relaxed within the time scale due to two- 
body collisions ((2.3(1 . This means that the stellar polytrope is no longer the equilibrium but a quasi- 
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/ \ plane. The thick arrows denote the evolution- 

P c' Pe' ary tracks in each simulation run (see Sec El. 



equilibrium state. In this section, we report the results of A^-body simulations |17j which were carried 
out to investigate how the stellar polytrope actually evolves. 

The iV-body experiment considered here is the same as that investigated in classic papers ([3101; 
see also That is, we confine the N particles interacting via Newtonian gravity in a spherical 

adiabatic wall, which reverses the radial components of the velocity if the particles reach the wall. 
Without loss of generality, we set the units as G = M = r e — 1. Note that the typical time scales 
appearing in this system are the dynamical time idyn <|2.2|1 and the global relaxation time driven by 
the two-body encounter, t xe \ J22J), which are basically scaled as td yn ~ 1 and t 10 \ ~ 0.1N/ \nN in 
our units. Table mi summarizes the five simulation runs. Here, in addition to the stellar polytropic 
initial state, we also consider the non-stellar polytropic state of the Hernquist model |37) . which was 
originally introduced to account for the empirical law of observed elliptical galaxies^. 

As was expected, the numerical simulations reveal that the stellar polytropic distribution gradually 
changes with time, on the time scale of two-body relaxation l|2.3|l . Furthermore, focusing on the 
evolutionary sequence, we found that the transient state starting from the initial stellar polytrope can 
be remarkably characterized by a sequence of stellar polytropes (run A, Bl, and B2). This is even true 
in the case starting from the Hernquist model (run CI). 

Let us show representative results taken from run A fFig l4.2(l . Figure l4~2T a~) plots snapshots of 
the density profile p(r), while Fig |4.2f b1 represents the distribution function f(e) as a function of the 
specific energy e = ^v 2 + <P(x). Note, for illustrative purpose, that each output result is artificially 
shifted to the two digits below. Only the final output with T = 400 represents the correct scales. 
In each figure, solid lines mean the initial stellar polytrope with n = 3 and the other lines indicate 



Table 4.1. Initial distributions and their evolutionary states 



run no. 


initial distribution 




parameters 


no. of particles 


transient state 


final state 


A 


stellar polytrope(n = 


3) 


D = 10, 000 


2,048 


stellar polytrope 


collapse 


Bl 


stellar polytrope(n = 


6) 


D = 110 


2,048 


stellar polytrope 


collapse 


B2 


stellar polytrope(n = 


6) 


D = 10, 000 


2,048 


stellar polytrope 


collapse 


CI 


Hernquist model 




a/r e = 0.5 


8,192 


stellar polytrope 


collapse 


C2 


Hernquist model 




a/r e = 1.0 


8,192 


none 


isothermal 
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the fitting results to the stellar polytrope by varying the polytrope index n. Note that the number 
of fitting parameters just dexreases to one, i.e. the polytrope index, since the total energy is well- 
conserved in the present situation. Fig 14. 21 shows that while the system gradually deviates from the 
initial polytropic state, the transient state still follows a sequence of stellar polytropes. The fitting 
results are remarkably good until the time exceeds T ~ 400, corresponding to 15i re i- After that, the 
system enters the gravothermally unstable regime and finally undergoes the core collapse. 

Now, focus on the evolutionary track in each simulation run summarized in the energy-radius-density 
contrast plane (Fig l4.1|l . where the filled circle represents the initial stellar polytrope. Interestingly, the 
density contrast of the transient state in run A initially decreases, but it eventually turns to increase. 
The turning point roughly corresponds to the stellar polytrope with index n ~ 5 — 6. Note, however, 
that the time evolution of the polytrope index itself is a monotonically increasing function of time |17j . 
This is indeed true for the other cases, indicating the Boltzmann iJ-theorem that any self-gravitating 
system tends to approach the isothermal state. A typical example is run C2, which finally reaches 
the stable isothermal state. However, as already shown in run A, not all the systems can reach the 
isothermal state. Fig l4. II indicates that no isothermal state is possible for a fixed value A > 0.335 |5P>] . 
which can be derived from the peak value of the trajectory. Further, stable stellar polytropes cease to 
exist at a high density contrast, D > D c rit- In fact, our simulations starting from the stellar polytropes 
finally underwent core collapse due to the gravothermal catastrophe (runs A, Bl, and B2). Though 
it might not be rigorously correct, the predicted value -D cr it provides a crude approximation to the 
boundary between stability and instability. 

Fig l4.3l plots snapshots of the distribution function taken from the other runs. The initial density 
contrast in run BI (Fig |4.3f a^ is relatively low (D — 110) and therefore the system slowly evolves 
by following a sequence of stellar polytropes. After T = 2000 ~ 74 t Te \, the system begins to deviate 
from the stable equilibrium sequence, leading to core collapse. Another noticeable case is run CI 
(Fig l4.3f b)). The Hcrnquist model as the initial distribution of run C has a cuspy density profile, 
p(r) oc l/r/(r + a) 3 , which behaves as p cx at the inner part |37) . The resultant distribution 
function f(e) shows singular behavior at the negative energy region, which cannot be described by 
the power-law distribution. Soon after, however, gravothermal expansion[3B] takes place and a flatter 
core is eventually formed. Then the system settles into a sequence of stellar polytropes and can be 
approximately described by the stellar polytrope with index n — 20 for a long time. Thus the stellar 
polytrope can be regarded as a quasi-attractor and a quasi-equilibrium state. 
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5 Conclusion and Discussion 

In this paper, we discussed issues arising from the Tsallis entropy for the thermodynamic properties of 
stellar self-gravitating systems, with a particular emphasis on the standard framework using normalized 
g-expectation values. It turns out that the new extremum-entropy state essentially remains unchanged 
from previous studies and is characterized by the stellar polytrope, although the distribution function 
shows several distinct properties. By considering these facts carefully, the thermodynamic temperature 
of the cxtrcmum state was identified through the modified Clausius relation and the specific heat was 
evaluated explicitly. A detailed analysis of the behavior of specific heat finally led to the conclusion 
that the onset of gravothermal instability remains unchanged with respect to choice of the statistical 
average for a system confined by an adiabatic wall (micro-canonical case) . As for a system surrounded 
by a thermal wall (canonical case), although the analysis has been skipped in this article, the stability 
of the system drastically depends on the choice of the statistical average |15Ulf>| . The existence of these 
thermodynamic instabilities can also be deduced rigorously from the variation of the entropy and free 
energy |14U15lll6j . As a result, above certain critical values of A or D, the thermodynamic instability 
appears at n > 5 for a system confined by an adiabatic wall. 

We performed a set of numerical simulations of long-term stellar dynamical evolution away from 
the isothermal state and found that the transient state of a system confined by an adiabatic wall can 
be remarkably fitted by a sequence of stellar polytropes. This is even true for the case in which the 
outer boundary in removed |17j . Therefore, the stellar polytropic distribution can be a quasi-attractor 
and a quasi-equilibrium state of a self-gravitating system. 
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A Appendix: Evaluation of the relaxation time 

Here we briefly evaluate the relaxation time for a self-gravitating many-body system l|2.3|l . For a more 
precise derivation, consult [1I18| . In the kinetic theory, a time scale of relaxation due to two-body 
collisions is given by 

trel ~ —, (A.l) 

anv 
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where n is the mean number density and v the average relative velocity of particles. In order to evaluate 
theamplitude of a two-body collision by the mutual gravitational force, we introduce the gravitational 
radius r g as follows, 

9 Gm 2 Gm , , . 

TTIV ~ > TV ~ —5-. (A. 2) 

If the impact parameter of the collision between two particles of identical mass m becomes smaller 
than the gravitational radius, i.e. b < r g , orbits of particles are significantly changed by the close 
encounter. It follows that the cross section of a two-body collision is given by 

a ~ Airr 2 . (A.3) 
By substituting this estimate and n ~ N/R 3 into (|A.1J) . we obtain 

v 3 R 3 

tlcl ~ Att G 2 m 2 N In (R/r g ) ' { ' 

where R and N are the system size and the number of particles, respectively. In the above expression 
(|A.4fl . we have included the so-called Coulomb logarithm term, \n(R/r g ), which appears due to the 
long-range nature of the interaction JTJ and is well-known in plasma theory. By taking its ratio to the 
dynamical time td yn ~ R/v, we obtain 

N 

tlcl ~4nh[N tdyn ' (A - 5) 
where the average velocity is estimated as v 2 ~ GNm/R from the virial theorem 12.4JI . 
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